
# REPLICATIONS CODE for “Did Border Closures Slow SARS-CoV-2?”

# Manuscript Authors: Mary A. Shiraef, Paul Friesen, Lukas Feddern, Mark A. Weiss, and COBAP Team
# Submited: December 12, 2021
# Code written by Paul Friesen (pfriesen@nd.edu) and Lukas Feddern (feddernlukas@gmail.com)
# Instructions - This file generates replication analysis and figures for the manuscript “Did Border Closures Slow SARS-CoV-2?”. The data files can be found in the same folder. Please make sure that the working directory is set to the replication folder. Please contact the code authors with questions or errors. 




# Set working Directory

setwd("~/X/X")

# Load Packages -----------------------

ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg))
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}


packages <- c("zoo",
              "rio",
              "tidyverse",
              "ggfortify",
              "stats",
              "astsa",
              "texreg",
              "forecast",
              "lubridate", 
              "urca",
              "reshape",
              "data.table",
              "readxl",
              "readr",
              "foreign",
              "tscount",
              "reshape2",
              "nortest",
              "car",
              "patchwork",
              "scales",
              "countrycode",
              "LambertW",
              "ggplot2",
              "plm",
              "PanelMatch")

ipak(packages)

# Descriptive Analysis over time -----------------
  # Data Management-----------------------------------

#Import data

COVID <- fread("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv")

COBAP_P <- fread("policy_list.csv")

##Time frame 

start_date <- lubridate::mdy("01/31/2020")
end_date <- lubridate::mdy("04/19/2021")
date <- seq(start_date, end_date, by="days")

##COBAP 

COBAP_P <- COBAP_P %>% 
  dplyr::select(ISO3, POLICY_TYPE:WORK_EXCEP) %>% 
  mutate(START_DATE = gsub("_", "-", START_DATE)) %>% 
  mutate(END_DATE = gsub("_", "-", END_DATE)) %>% 
  mutate(START_DATE = mdy(START_DATE)) %>% 
  mutate(END_DATE = mdy(END_DATE)) %>% 
  mutate(DURATION = END_DATE-START_DATE)


#Get country list

COVID <- COVID %>% 
  pivot_longer(!(`Province/State`|`Country/Region`|Lat|Long), names_to = "date", values_to = "cases") %>%
  mutate(date = mdy(date)) %>% 
  mutate(new_cases = cases - shift(cases, 1L, type="lag")) %>% 
  mutate(ISO3 = countrycode(`Province/State`, origin = 'country.name', destination = 'iso3c')) %>% 
  mutate(ISO3 = coalesce(ISO3, countrycode(`Country/Region`, origin = 'country.name', destination = 'iso3c')))

COVID$ISO3[COVID$`Country/Region`=="Kosovo"] <- "XK"

country_list <- COVID %>%   
  dplyr::select(ISO3, cases) %>% 
  drop_na() %>% 
  dplyr::select(ISO3) %>% 
  unique()


COVID <- COVID %>%   
  dplyr::select(ISO3, date, new_cases) %>% 
  group_by(ISO3, date) %>% 
  summarise(
    new_cases = sum(new_cases, na.rm = T)) %>% 
  ungroup() %>% 
  filter(between(date, start_date, end_date))

#Create time series data set

COBAP.daily <- as_tibble(country_list) %>% 
  slice(rep(1:n(), each = length(date))) %>% 
  cbind(date) %>% 
  left_join(COBAP_P, by = c("ISO3"="ISO3", "date"="START_DATE")) %>% 
  mutate(complete_new = if_else(POLICY_TYPE == "COMPLETE", 1, 0)) %>% 
  mutate(partial_new = if_else(POLICY_TYPE == "PARTIAL", 1, 0)) %>%
  mutate(borderclosure_new = if_else(POLICY_SUBTYPE == "BORDER_CLOSURE", 1, 0)) %>% 
  mutate(partial2_new = if_else(POLICY_TYPE == "PARTIAL" & POLICY_SUBTYPE != "VISA_BAN", 1, 0))


COBAP.daily$partial_new = dplyr::recode(COBAP.daily$partial_new, .missing = 0)
COBAP.daily$partial2_new = dplyr::recode(COBAP.daily$partial2_new, .missing = 0)
COBAP.daily$complete_new = dplyr::recode(COBAP.daily$complete_new, .missing = 0)

#Count policies

COBAP.daily <- COBAP.daily %>% 
  group_by(ISO3) %>% 
  mutate(partial = purrr::accumulate(partial_new, `+`)) %>% 
  mutate(partial2 = purrr:: accumulate(partial2_new, `+`)) %>% 
  mutate(complete = purrr::accumulate(complete_new, `+`)) %>% 
  ungroup() %>% 
  mutate(air_c = 0) 

#Count policies for each country recorded

COBAPdata_list = split(COBAP.daily, COBAP.daily$ISO3)

for (i in seq_along(COBAPdata_list)) {
  for (j in 1:nrow(COBAPdata_list[[i]])) {
    if (COBAPdata_list[[i]][j,33] ==1 && is.na(COBAPdata_list[[i]][j,5])==F){
      COBAPdata_list[[i]]$partial2[COBAPdata_list[[i]]$date >= COBAPdata_list[[i]][j,5]$END_DATE] <- COBAPdata_list[[i]]$partial2[COBAPdata_list[[i]]$date >= COBAPdata_list[[i]][j,5]$END_DATE]-1
    }else if (COBAPdata_list[[i]][j,30] ==1 && is.na(COBAPdata_list[[i]][j,5])==F){
      COBAPdata_list[[i]]$complete[COBAPdata_list[[i]]$date >= COBAPdata_list[[i]][j,5]$END_DATE] <- COBAPdata_list[[i]]$complete[COBAPdata_list[[i]]$date >= COBAPdata_list[[i]][j,5]$END_DATE]-1
    }else{
      COBAPdata_list[[i]]$complete <- COBAPdata_list[[i]]$complete  
    }
  }
}

#Summarize count for each day 

for (i in seq_along(COBAPdata_list)){
  COBAPdata_list[[i]]<-COBAPdata_list[[i]] %>% 
    group_by(date) %>% 
    summarise(
      partial_closure = max(partial2),
      complete_closure = max(complete),
      ISO3 = ISO3) %>% 
    ungroup() %>% 
    distinct()
}


#Create global dataset
COBAP.daily <- merge_all(COBAPdata_list)

COVID <- COVID %>%   
  dplyr::select(ISO3, date, new_cases) %>% 
  drop_na(ISO3)

#Merge case dataset and add policy count
COBAP.daily.global <- COBAP.daily %>% 
  dplyr::select(ISO3, date, complete_closure, partial_closure) %>% 
  left_join(COVID, by = c("date"="date", "ISO3"="ISO3")) %>% 
  group_by(date) %>% 
  summarise(
    global_complete = sum(complete_closure, na.rm = T),
    global_borderclosures = sum(partial_closure, na.rm = T),
    global_new_cases = sum(new_cases, na.rm = T)
  ) 

COBAP.daily.global <- COBAP.daily.global %>% 
  mutate(day = wday(COBAP.daily.global$date)) %>% 
  mutate(week = week(COBAP.daily.global$date))

#Correction of outlier (Dec-10-2020: Source https://web.archive.org/web/20210618230847/https://ourworldindata.org/grapher/daily-covid-cases-deaths?country=~OWID_WRL)

COBAP.daily.global$global_new_cases<-replace(COBAP.daily.global$global_new_cases, 315, 674901)

  # Data Analysis and Plotting ----------------------------

# Plots
plot_global.cases <- COBAP.daily.global %>% 
  ggplot()+
  geom_line(mapping = aes(x = date, y = global_new_cases), color = "red") +
  ggtitle("Global Daily New COVID-19 Cases") +
  ylab("New COVID-19 Cases")+
  xlab("Date")+
  theme_bw()

plot_global.complete <- COBAP.daily.global %>% 
  ggplot()+
  geom_line(mapping = aes(x = date, y = global_complete), color = "blue")+
  ggtitle("Global Number of Complete Closures") +
  ylab("Number of Complete Closures")+
  xlab("Date")+
  theme_bw()

vizdf <- COBAP.daily.global %>% 
  dplyr::select(date, global_borderclosures, global_complete)

vizdf <- melt(setDT(vizdf), id="date")


# Create Figure 1
ggplot() + 
  geom_area(data = vizdf, aes(x=date, y=value, fill=variable))+
  geom_line(data = COBAP.daily.global, mapping = aes(x = date, y = global_new_cases/1438.793), color = "black")+
  scale_y_continuous(
    name ="Total Number of Policies (Border closures)",
    sec.axis = sec_axis(~.*1438.793,name="New COVID-19 Cases", labels = comma)
  )+
  theme_bw()+
  scale_fill_discrete(name = "Border closure type", labels = c("Partial border closure \n (excl. visa bans)", "Complete closure"))+
  labs(x = "Date") +
  scale_fill_manual(
    name = "Border closure type",
    values = c("#A9A9A9", "darkred"),
    labels = c("Partial border closure \n (excl. visa bans)", "Complete closure")
  )



# Panel Matching Set-Up -------------------------------
  # Load Datasets --------------

COBAP <-  read_excel("COBAP_timeseries_Dec_11_2021.xlsx")

COBAP_Island_Countries_List <- read_excel("Islands_list.xlsx", col_names = F)

ECON <- read_excel("data_with_sero_ifr_econ_pol.xlsx")

JH <- read_csv("JH_Cases.csv")

CGRT_School <- read_excel("OxCGRT_timeseries_all.xlsx",  sheet = "c1_school_closing")
  
CGRT_Workplace <- read_excel("OxCGRT_timeseries_all.xlsx",  sheet = "c2_workplace_closing")

CGRT_PublicEvents <- read_excel("OxCGRT_timeseries_all.xlsx",  sheet = "c3_cancel_public_events")

CGRT_Gatherings <- read_excel("OxCGRT_timeseries_all.xlsx",  sheet = "c4_restrictions_on_gatherings")

CGRT_Transport <- read_excel("OxCGRT_timeseries_all.xlsx",  sheet = "c5_close_public_transport")

CGRT_StayHome <- read_excel("OxCGRT_timeseries_all.xlsx",  sheet = "c6_stay_at_home_requirements")

CGRT_Movement <- read_excel("OxCGRT_timeseries_all.xlsx",  sheet = "c7_movementrestrictions")



  # Johns Hopkins Data Cleaning ----------------------
JH <- JH[ , -c(1,3,4)]
JH_country <-  aggregate(.~`Country/Region`, JH, sum)

JH_country$iso3c <- c("AFG", "ALB", "DZA", "AND", "AGO", "ATG", "ARG", "ARM", "AUS", "AUT", "AZE", "BHS", "BHR", "BGD", "BRB", "BLR", "BEL", "BLZ", "BEN", "BTN", "BOL", "BIH", "BWA", "BRA", "BRN", "BGR", "BFA", "MMR", "BDI", "CPV", "KHM", "CMR", "CAN", "CAF", "TCD", "CHL", "CHN", "COL", "COM", "COG", "COD", "CRI", "CIV", "HRV", "CUB", "CYP", "CZE", "DNK", "0", "DJI", "DMA", "DOM", "ECU", "EGY", "SLV", "GNQ", "ERI", "EST", "SWZ", "ETH", "FJI", "FIN", "FRA", "GAB", "GMB", "GEO", "DEU", "GHA", "GRC", "GRD", "GTM", "GIN", "GNB", "GUY", "HTI", "VAT", "HND", "HUN", "ISL", "IND", "IDN", "IRN", "IRQ", "IRL", "ISR", "ITA", "JAM", "JPN", "JOR", "KAZ", "KEN", "KIR", "KOR", "XKX", "KWT", "KGZ", "LAO", "LVA", "LBN", "LSO", "LBR", "LBY", "LIE", "LTU", "LUX", "MDG", "MWI", "MYS", "MDV", "MLI", "MLT", "MHL", "MRT", "MUS", "MEX", "FSM", "MDA", "MCO", "MNG", "MNE", "MAR", "MOZ", "0", "NAM", "NPL", "NLD", "NZL", "NIC", "NER", "NGA", "MKD", "NOR", "OMN", "PAK", "PLW", "PAN", "PNG", "PRY", "PER", "PHL", "POL", "PRT", "QAT", "ROU", "RUS", "RWA", "KNA", "LCA", "VCT", "WSM", "SMR", "STP", "SAU", "SEN", "SRB", "SYC", "SLE", "SGP", "SVK", "SVN", "SLB", "SOM", "ZAF", "SSD", "ESP", "LKA", "SDN", "0", "SUR", "SWE", "CHE", "SYR", "TWN", "TJK", "TZA", "THA", "TLS", "TGO", "TNO", "TTO", "TUN", "TUR", "UGA", "UKR", "ARE", "GBR", "URY", "USA", "UZB", "VUT", "VEN", "VNM", "PSE", "YEM", "ZMB", "ZWE")

JH_long <- reshape::melt(JH_country, id.vars=c("iso3c", "Country/Region"))

#June5
JH_long <- JH_long[1:98196, ]

JH_long$time_count <- rep(22:522,each=196)

JH_long$week_time_count <- JH_long$time_count-4

JH_long$week_count <- ifelse(JH_long$week_time_count>0, round((JH_long$week_time_count/7)+0.499,0), NA)

JH_long <- JH_long[-which(JH_long$iso3c=="0"),]
JH_long$iso3c <- as.character(JH_long$iso3c)

JH_Cleaned <- aggregate(value ~ iso3c + week_count, JH_long, max)

rm(JH, JH_country, JH_long)

  # Economist Data Cleaning --------------------------------

ECON <- ECON[which(ECON$`C/T`=="c"), ]

ECON <- ECON[-c(28517:29128), ]

ECON$year <-as.numeric(substr(ECON$date, start=1, stop=4))
ECON$month <- as.numeric(substr(ECON$date, start=6, stop=7))
ECON$day <- as.numeric(substr(ECON$date, start=9, stop=10))

ECON$time_count <- ((ECON$year-2020)*366) + (ifelse(ECON$month==1, 0, ifelse(ECON$month==2, 31, ifelse(ECON$month==3, 60, ifelse(ECON$month==4, 91, ifelse(ECON$month==5, 121, ifelse(ECON$month==6, 152, ifelse(ECON$month==7, 182, ifelse(ECON$month==8, 213, ifelse(ECON$month==9, 244, ifelse(ECON$month==10, 274, ifelse(ECON$month==11, 305, 335)))))))))))) + ECON$day

ECON$week_time_count <- ECON$time_count-4

ECON$week_count <- ifelse(ECON$week_time_count>0, round((ECON$week_time_count/7)+0.499,0), NA)

ECON_E <- ECON[,c(2,4,6,9:15,22,124,129:136,143,146,147,158:163)]

ECON_E$gdp_per_capita <- as.numeric(ECON_E$gdp_per_capita)
ECON_E$human_development_index <- as.numeric(ECON_E$human_development_index)
ECON_E$v2x_libdem <- as.numeric(ECON_E$v2x_libdem)
ECON_E$v2x_freexp_altinf <- as.numeric(ECON_E$v2x_freexp_altinf)
ECON_E$median_age <- as.numeric(ECON_E$median_age)
ECON_E$aged_70_older <- as.numeric(ECON_E$aged_70_older)
ECON_E$life_expectancy <- as.numeric(ECON_E$life_expectancy)
ECON_E$daily_tests_per_100k <- as.numeric(ECON_E$daily_tests_per_100k)
ECON_E$area_ifr <- as.numeric(ECON_E$area_ifr)
ECON_E$pop.2020 <- as.numeric(ECON_E$pop.2020)
ECON_E$first.dose.pct <- as.numeric(ECON_E$first.dose.pct)
ECON_E$second.dose.pct <- as.numeric(ECON_E$second.dose.pct)


# Combining variables into week aggregations - different variable transformation operation

ECON_E_fixed <- ECON_E[,c(1:10,12:13,29)]

ECON_E_sum <- ECON_E[,c(1,11,22:23,29)]

ECON_E_fixed <- ECON_E_fixed[!duplicated(ECON_E_fixed),]

ECON_E_sum <- ECON_E_sum %>%
  group_by(iso3c, week_count) %>% 
  summarise_each(funs(sum))  

# Recombine week aggregations
ECON_Cleaned <- merge(ECON_E_sum, ECON_E_fixed, by = c("iso3c", "week_count"))

rm(ECON, ECON_E, ECON_E_fixed, ECON_E_sum)



  # CGR Data Cleaning  --------------------

CGRT_School <- CGRT_School[, -2]
CGRT_Workplace <- CGRT_Workplace[, -2]
CGRT_PublicEvents <- CGRT_PublicEvents[, -2]
CGRT_Gatherings <- CGRT_Gatherings[, -2]
CGRT_Transport <- CGRT_Transport[, -2]
CGRT_StayHome <- CGRT_StayHome[, -2]
CGRT_Movement <- CGRT_Movement[, -2]

# Convert variables to long form
CGRT_School <- reshape2::melt(CGRT_School, id.vars=c("country_code"))
CGRT_Workplace <- reshape2::melt(CGRT_Workplace, id.vars=c("country_code"))
CGRT_PublicEvents <- reshape2::melt(CGRT_PublicEvents, id.vars=c("country_code"))
CGRT_Gatherings <- reshape2::melt(CGRT_Gatherings, id.vars=c("country_code"))
CGRT_Transport <- reshape2::melt(CGRT_Transport, id.vars=c("country_code"))
CGRT_StayHome <- reshape2::melt(CGRT_StayHome, id.vars=c("country_code"))
CGRT_Movement <- reshape2::melt(CGRT_Movement, id.vars=c("country_code"))

# Combine into single dataset
CGRT <- cbind(CGRT_School, CGRT_Workplace,  CGRT_PublicEvents, CGRT_Gatherings, CGRT_Transport, CGRT_StayHome, CGRT_Movement)

# Remove redundant columns
CGRT <- CGRT[ , -c(4,5,7,8,10,11,13,14,16,17,19,20)]
colnames(CGRT) <- c("country", "date", "school", "workplace", "publicevent", "gatherings", "transport", "stayhome", "movement")

# Create numeric date columns
CGRT$day <- as.numeric(substr(CGRT$date, start=1, stop=2))
CGRT$month <- substr(CGRT$date, start=3, stop=5)
CGRT$year <- as.numeric(substr(CGRT$date, start=6, stop=9))

CGRT$month <- ifelse(CGRT$month=="Jan", 1, ifelse(CGRT$month=="Feb", 2, ifelse(CGRT$month=="Mar", 3, ifelse(CGRT$month=="Apr", 4, ifelse(CGRT$month=="May", 5, ifelse(CGRT$month=="Jun", 6, ifelse(CGRT$month=="Jul", 7, ifelse(CGRT$month=="Aug", 8, ifelse(CGRT$month=="Sep", 9, ifelse(CGRT$month=="Oct", 10, ifelse(CGRT$month=="Nov", 11, ifelse(CGRT$month=="Dec", 12, NA))))))))))))

# Create day count variables
CGRT$time_count <- ((CGRT$year-2020)*366) + (ifelse(CGRT$month==1, 0, ifelse(CGRT$month==2, 31, ifelse(CGRT$month==3, 60, ifelse(CGRT$month==4, 91, ifelse(CGRT$month==5, 121, ifelse(CGRT$month==6, 152, ifelse(CGRT$month==7, 182, ifelse(CGRT$month==8, 213, ifelse(CGRT$month==9, 244, ifelse(CGRT$month==10, 274, ifelse(CGRT$month==11, 305, 335)))))))))))) + CGRT$day

# Adjust for start date
CGRT$week_time_count <- CGRT$time_count-4

# Create week counts
CGRT$week_count <- ifelse(CGRT$week_time_count>0, round((CGRT$week_time_count/7)+0.499,0), NA)

# Standardized variables between 0 and 1
CGRT$school_i <- CGRT$school/3
CGRT$workplace_i <- CGRT$workplace/3
CGRT$publicevent_i <- CGRT$publicevent/2
CGRT$gatherings_i <- CGRT$gatherings/4
CGRT$transport_i <- CGRT$transport/2
CGRT$stayhome_i <- CGRT$stayhome/3
CGRT$movement_i <- CGRT$movement/2

# Combine indicators into new Stringency Index
CGRT$Stringency_revised <- (CGRT$school_i + CGRT$workplace_i + CGRT$publicevent_i + CGRT$gatherings_i + CGRT$transport_i + CGRT$stayhome_i + CGRT$movement_i)/7


# Remove redundant observations
CGRT_E <- CGRT[745:132246, c(1,15:23)]

# Create weekly average statistics
CGRT_Cleaned <- CGRT_E %>%
  group_by(country, week_count) %>% 
  summarise_each(funs(mean))  


rm(CGRT, CGRT_E, CGRT_Gatherings, CGRT_Movement, CGRT_PublicEvents, CGRT_School, CGRT_StayHome, CGRT_Transport, CGRT_Workplace)


  # Combining Datasets --------------

colnames(CGRT_Cleaned)[1] <- "iso3c"
colnames(COBAP)[1] <- "iso3c"

COBAP$week_count <- rep(4:67, 252)

DF_Full <- merge(COBAP, ECON_Cleaned, by=c("iso3c", "week_count"), all=T)

DF_Full <- merge(DF_Full, CGRT_Cleaned[, c(1,2,10)], by=c("iso3c", "week_count"), all=T)

DF_Full <- merge(DF_Full, JH_Cleaned, by=c("iso3c", "week_count"), all=T)

DF_Full <- DF_Full[which(is.na(DF_Full$week_count)==F), ]
DF_Full <- DF_Full[-which(DF_Full$week_count>=68), ]


  # Variable Transformation ----------------
DF_Full$completenew <- as.numeric(DF_Full$completenew)
DF_Full$completenew_D <- ifelse(DF_Full$completenew>=1, 1, 0)
DF_Full$completenew_D <- ifelse(DF_Full$week_count<=3, 0, DF_Full$completenew_D)

DF_Full$partialnew <- as.numeric(DF_Full$partialnew)
DF_Full$partialnew_D <- ifelse(DF_Full$partialnew>=1, 1, 0)
DF_Full$partialnew_D <- ifelse(DF_Full$week_count<=3, 0, DF_Full$partialnew_D)

DF_Full$specific_countrynew <- as.numeric(DF_Full$specific_countrynew)
DF_Full$specific_countrynew_D <- ifelse(DF_Full$specific_countrynew>=1, 1, 0)
DF_Full$specific_countrynew_D <- ifelse(DF_Full$week_count<=3, 0, DF_Full$specific_countrynew_D)

DF_Full$work_excepnew <- as.numeric(DF_Full$work_excepnew)
DF_Full$work_excepnew_D <- ifelse(DF_Full$work_excepnew>=1, 1, 0)
DF_Full$work_excepnew_D <- ifelse(DF_Full$week_count<=3, 0, DF_Full$work_excepnew_D)

DF_Full$citizen_excepnew <- as.numeric(DF_Full$citizen_excepnew)
DF_Full$citizen_excepnew_D <- ifelse(DF_Full$citizen_excepnew>=1, 1, 0)
DF_Full$citizen_excepnew_D <- ifelse(DF_Full$week_count<=3, 0, DF_Full$citizen_excepnew_D)

DF_Full$essential_onlynew <- as.numeric(DF_Full$essential_onlynew)
DF_Full$essential_onlynew_D <- ifelse(DF_Full$essential_onlynew>=1, 1, 0)
DF_Full$essential_onlynew_D <- ifelse(DF_Full$week_count<=3, 0, DF_Full$essential_onlynew_D)

DF_Full$visa_bannew <- as.numeric(DF_Full$visa_bannew)
DF_Full$visa_bannew_D <- ifelse(DF_Full$visa_bannew>=1, 1, 0)
DF_Full$visa_bannew_D <- ifelse(DF_Full$week_count<=3, 0, DF_Full$visa_bannew_D)

DF_Full$citizenship_bannew <- as.numeric(DF_Full$citizenship_bannew)
DF_Full$citizenship_bannew_D <- ifelse(DF_Full$citizenship_bannew>=1, 1, 0)
DF_Full$citizenship_bannew_D <- ifelse(DF_Full$week_count<=3, 0, DF_Full$citizenship_bannew_D)

DF_Full$history_bannew <- as.numeric(DF_Full$history_bannew)
DF_Full$history_bannew_D <- ifelse(DF_Full$history_bannew>=1, 1, 0)
DF_Full$history_bannew_D <- ifelse(DF_Full$week_count<=3, 0, DF_Full$history_bannew_D)

DF_Full$border_closurenew <- as.numeric(DF_Full$border_closurenew)
DF_Full$border_closurenew_D <- ifelse(DF_Full$border_closurenew>=1, 1, 0)
DF_Full$border_closurenew_D <- ifelse(DF_Full$week_count<=3, 0, DF_Full$border_closurenew_D)

DF_Full <- DF_Full %>%
  group_by(iso3c) %>%
  fill(human_development_index) %>%
  fill(human_development_index, .direction = "up")

DF_Full <- DF_Full %>%
  group_by(iso3c) %>%
  fill(v2x_libdem) %>%
  fill(v2x_libdem, .direction = "up")

DF_Full <- DF_Full %>%
  group_by(iso3c) %>%
  fill(v2x_freexp_altinf) %>%
  fill(v2x_freexp_altinf, .direction = "up")

DF_Full <- DF_Full %>%
  group_by(iso3c) %>%
  fill(pop.2020) %>%
  fill(pop.2020, .direction = "up")

DF_Full <- DF_Full %>%
  group_by(iso3c) %>%
  fill(aged_70_older) %>%
  fill(aged_70_older, .direction = "up")

DF_Full <- DF_Full %>%
  group_by(iso3c) %>%
  fill(life_expectancy) %>%
  fill(life_expectancy, .direction = "up")

DF_Full <- DF_Full %>%
  group_by(iso3c) %>%
  fill(gdp_per_capita) %>%
  fill(gdp_per_capita, .direction = "up")

DF_Full <- DF_Full %>%
  group_by(iso3c) %>%
  fill(area_ifr) %>%
  fill(area_ifr, .direction = "up")

DF_Full <- DF_Full %>%
  group_by(iso3c) %>%
  fill(median_age) %>%
  fill(median_age, .direction = "up")


nocases_early <- DF_Full[which(DF_Full$week_count==3 & DF_Full$value==0), ]$iso3c
zero2 <- intersect(which(DF_Full$week_count==2), which(DF_Full$iso3c %in% nocases_early==TRUE))
zero1 <- intersect(which(DF_Full$week_count==1), which(DF_Full$iso3c %in% nocases_early==TRUE))

DF_Full[zero1, which(colnames(DF_Full)=="value")] <- 0    
DF_Full[zero2, which(colnames(DF_Full)=="value")] <- 0            
DF_Full <- DF_Full[-which(duplicated(DF_Full[,1:2])==T), ]

DF_Full <- pdata.frame(DF_Full, index = c("iso3c", "week_count"))

DF_Full$value_1 <- plm::lag(DF_Full$value,1)
DF_Full$new_cases_UP <- DF_Full$value - DF_Full$value_1

DF_Full$new_cases_UP_100k <-  DF_Full$new_cases_UP/(DF_Full$pop.2020/100000)

DF_Full$new_cases_UP_100k <- ifelse(DF_Full$new_cases_UP_100k<0, 0, DF_Full$new_cases_UP_100k)

DF_Full$new_cases_UP_100k_1 <- plm::lag(DF_Full$new_cases_UP_100k,1)
DF_Full$new_cases_UP_100k_change <- DF_Full$new_cases_UP_100k-DF_Full$new_cases_UP_100k_1
DF_Full <- DF_Full[-which(is.na(DF_Full$new_cases_UP_100k_change)==T), ]


DF_Full$Stringency_revised_1 <- plm::lag(DF_Full$Stringency_revised, 1)
DF_Full$stringency_revised_change <- DF_Full$Stringency_revised - DF_Full$Stringency_revised_1
DF_Full$Stringency_Treat <- ifelse(DF_Full$stringency_revised_change>=0.25, 1, 0)

DF_Full <- as.data.frame(DF_Full)
DF_Full$countryN <- as.numeric(as.factor(DF_Full$iso3c))
DF_Full$week_count <- as.integer(DF_Full$week_count)

DF_Full$pop.2020_LOG <- log(DF_Full$pop.2020)
DF_Full$area_ifr_LOG <- log(DF_Full$area_ifr)
DF_Full$gdp_per_capita_LOG <- log(DF_Full$gdp_per_capita)

DF_Full[which(is.na(DF_Full$first.dose.pct)==T),]$first.dose.pct <- 0
DF_Full[which(is.na(DF_Full$second.dose.pct)==T),]$second.dose.pct <- 0

ihs <- function(x) {
  y <- log(x + sqrt(x^2 + 1))
  return(y)
}

DF_Full$new_cases_UP_100k_change_ihs <- ihs(DF_Full$new_cases_UP_100k_change)



  # Dependent Variable Diagnostics ----------------------

plot.new()

LambertW::test_norm(DF_Full$new_cases_UP_100k)

LambertW::test_norm(DF_Full$new_cases_UP_100k_change)

LambertW::test_norm(DF_Full$new_cases_UP_100k_change_ihs)

# Panel Match Analysis -----------------



  # Lockdowns ----------------------------

Lockdown_NN <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", refinement.method = "none", treatment = "Stringency_Treat", matching = TRUE, data = DF_Full, match.missing = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

pt_DL_NN <- placeboTest(Lockdown_NN, data.in = DF_Full, plot = F)
pt_DL_NN_t_3_est <- pt_DL_NN$estimates[1]
pt_DL_NN_t_2_est <- pt_DL_NN$estimates[2]

pt_DL_NN_t_3_high <- quantile(pt_DL_NN$bootstrapped.estimates[ ,1],0.975)
pt_DL_NN_t_2_high <- quantile(pt_DL_NN$bootstrapped.estimates[ ,2],0.975)
pt_DL_NN_t_3_low <- quantile(pt_DL_NN$bootstrapped.estimates[ ,1],0.025)
pt_DL_NN_t_2_low <- quantile(pt_DL_NN$bootstrapped.estimates[ ,2],0.025)

Results_DL_NN <- PanelEstimate(sets = Lockdown_NN, data = DF_Full)
summary(Results_DL_NN)


# Delete Stringency Revised from lag period
Lockdown_cbps <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "Stringency_Treat", refinement.method = "CBPS.match", covs.formula = ~ human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older +  median_age + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + first.dose.pct + second.dose.pct + daily_tests_per_100k + I(lag(daily_tests_per_100k, 1:3)) + Stringency_revised + I(lag(completenew_D, 1:3)) + completenew_D + I(lag(partialnew_D, 1:3)) + partialnew_D, data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

pt_DL_cbps <- placeboTest(Lockdown_cbps, data.in = DF_Full, plot = F)

pt_DL_cbps_t_3_est <- pt_DL_cbps$estimates[1]
pt_DL_cbps_t_2_est <- pt_DL_cbps$estimates[2]

pt_DL_cbps_t_3_high <- quantile(pt_DL_cbps$bootstrapped.estimates[ ,1],0.975)
pt_DL_cbps_t_2_high <- quantile(pt_DL_cbps$bootstrapped.estimates[ ,2],0.975)
pt_DL_cbps_t_3_low <- quantile(pt_DL_cbps$bootstrapped.estimates[ ,1],0.025)
pt_DL_cbps_t_2_low <- quantile(pt_DL_cbps$bootstrapped.estimates[ ,2],0.025)

Results_DL_cbps <- PanelEstimate(sets = Lockdown_cbps, data = DF_Full)
summary(Results_DL_cbps)







# Generate Figure 2

plot_table_NN <- rbind(pt_DL_NN_t_3_est, pt_DL_NN_t_2_est, 0, as.data.frame(Results_DL_NN$estimates))
colnames(plot_table_NN) <- "estimate"
plot_table_NN$high <- c(pt_DL_NN_t_3_high, pt_DL_NN_t_2_high, 0, Results_DL_NN$estimates + (Results_DL_NN$standard.error*1.96))
plot_table_NN$low <- c(pt_DL_NN_t_3_low, pt_DL_NN_t_2_low, 0, Results_DL_NN$estimates - (Results_DL_NN$standard.error*1.96))
plot_table_NN$Matching <- "No Matching" 
rownames(plot_table_NN)[1:4] <- c("t-3","t-2","t-1","t")
plot_table_NN$Week <- rownames(plot_table_NN)


plot_table_R <- rbind(pt_DL_cbps_t_3_est, pt_DL_cbps_t_2_est, 0, as.data.frame(Results_DL_cbps$estimates))
colnames(plot_table_R) <- "estimate"
plot_table_R$high <- c(pt_DL_cbps_t_3_high, pt_DL_cbps_t_2_high, 0, Results_DL_cbps$estimates + (Results_DL_cbps$standard.error*1.96))
plot_table_R$low <- c(pt_DL_cbps_t_3_low, pt_DL_cbps_t_2_low, 0, Results_DL_cbps$estimates - (Results_DL_cbps$standard.error*1.96))
plot_table_R$Matching <- "Matching + Refinement" 
rownames(plot_table_R)[1:4] <- c("t-3","t-2","t-1","t")
plot_table_R$Week <- rownames(plot_table_R)

plot_table_DL <- rbind(plot_table_R, plot_table_NN)
plot_table_DL$Week <- factor(plot_table_DL$Week, levels = c("t-3", "t-2", "t-1", "t", "t+1", "t+2", "t+3", "t+4", "t+5"))

plot_table_DL$`Estimation Strategy` <- factor(plot_table_DL$Matching, levels = c("No Matching", "Matching + Refinement"))



# Plot
pd = position_dodge(.6)   
ggplot(plot_table_DL, aes(x = Week, y = estimate, color = `Estimation Strategy`)) + geom_hline(yintercept = 0, color = "gray40", linetype="dashed") +
  geom_point(shape = 16, size  = 6, position = pd) +
  geom_errorbar(aes(ymin  = low, ymax  = high), width = 0, size  = 1.45, position = pd) +
  theme_bw() + theme(axis.title = element_text(face = "bold")) + ylab("Rate of Change - New Weekly Cases Per Capita (IHST)") + scale_color_manual(values= c("gray70", "#689950")) + theme(legend.text=element_text(size=20), axis.text = element_text(size = 20), axis.title = element_text(size = 20), text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + ggtitle("Domestic Lockdowns") + scale_y_continuous(limits = c(-1.5, 1.5), breaks=c(-1.5,-1,-0.5,0,0.5,1,1.5))






# Analysis of Refinement strategies - Change refinement method and re-run as needed. 

Lockdown_A <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "Stringency_Treat", matching = TRUE, refinement.method = "CBPS.match", covs.formula = ~ human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older +  median_age + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + first.dose.pct + second.dose.pct + daily_tests_per_100k + I(lag(daily_tests_per_100k, 1:3)) + Stringency_revised + I(lag(partialnew_D, 1:3)) + partialnew_D + I(lag(completenew_D, 1:3)) + completenew_D, data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = F)

cvb_LD_A <- get_covariate_balance(Lockdown_A$att,DF_Full, covariates = c("gdp_per_capita_LOG", "new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "median_age", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy", "partialnew_D", "Stringency_Treat", "completenew_D"))

Var <- ncol(cvb_LD_A)
cvb_LD_A_totals <- data.frame(sum(ifelse(abs(cvb_LD_A[1,])>=0.20, abs(cvb_LD_A[1,])-0.20, 0))/Var, sum(ifelse(abs(cvb_LD_A[2,])>=0.20, abs(cvb_LD_A[2,])-0.20, 0))/Var, sum(ifelse(abs(cvb_LD_A[3,])>=0.20, abs(cvb_LD_A[3,])-0.20, 0))/Var, sum(ifelse(abs(cvb_LD_A[4,])>=0.20, abs(cvb_LD_A[4,])-0.20, 0))/Var)
cvb_LD_A_totals


# Imbalance scores by refinement strategy
  # None - 0.0883
  # Mahal - 0.120
  # PS - 0.053
  # CBPS - 0.050

  # Complete Closures ------------

Complete_NN <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", refinement.method = "none", treatment = "completenew_D", matching = TRUE, data = DF_Full, match.missing = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

pt_CC_NN <- placeboTest(Complete_NN, data.in = DF_Full, plot = F)
pt_CC_NN_t_3_est <- pt_CC_NN$estimates[1]
pt_CC_NN_t_2_est <- pt_CC_NN$estimates[2]

pt_CC_NN_t_3_high <- quantile(pt_CC_NN$bootstrapped.estimates[ ,1],0.975)
pt_CC_NN_t_2_high <- quantile(pt_CC_NN$bootstrapped.estimates[ ,2],0.975)
pt_CC_NN_t_3_low <- quantile(pt_CC_NN$bootstrapped.estimates[ ,1],0.025)
pt_CC_NN_t_2_low <- quantile(pt_CC_NN$bootstrapped.estimates[ ,2],0.025)

Results_CC_NN <- PanelEstimate(sets = Complete_NN, data = DF_Full)

Complete_ps <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "completenew_D", refinement.method = "ps.match", covs.formula = ~ human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older +  median_age + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + first.dose.pct + second.dose.pct + daily_tests_per_100k + I(lag(daily_tests_per_100k, 1:3))+ I(lag(Stringency_revised, 1:3)) + Stringency_revised + I(lag(partialnew_D, 1:3)) + partialnew_D + I(lag(Stringency_Treat, 1:3)) + Stringency_Treat  , data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

pt_CC_ps <- placeboTest(Complete_ps, data.in = DF_Full, plot = F)

pt_CC_ps_t_3_est <- pt_CC_ps$estimates[1]
pt_CC_ps_t_2_est <- pt_CC_ps$estimates[2]

pt_CC_ps_t_3_high <- quantile(pt_CC_ps$bootstrapped.estimates[ ,1],0.975)
pt_CC_ps_t_2_high <- quantile(pt_CC_ps$bootstrapped.estimates[ ,2],0.975)
pt_CC_ps_t_3_low <- quantile(pt_CC_ps$bootstrapped.estimates[ ,1],0.025)
pt_CC_ps_t_2_low <- quantile(pt_CC_ps$bootstrapped.estimates[ ,2],0.025)

Results_CC_ps <- PanelEstimate(sets = Complete_ps, data = DF_Full)

summary(Results_CC_ps)



# Generate Figure 3

plot_table_NN <- rbind(pt_CC_NN_t_3_est, pt_CC_NN_t_2_est, 0, as.data.frame(Results_CC_NN$estimates))
colnames(plot_table_NN) <- "estimate"
plot_table_NN$high <- c(pt_CC_NN_t_3_high, pt_CC_NN_t_2_high, 0, Results_CC_NN$estimates + (Results_CC_NN$standard.error*1.96))
plot_table_NN$low <- c(pt_CC_NN_t_3_low, pt_CC_NN_t_2_low, 0, Results_CC_NN$estimates - (Results_CC_NN$standard.error*1.96))
plot_table_NN$Matching <- "No Matching" 
rownames(plot_table_NN)[1:4] <- c("t-3","t-2","t-1","t")
plot_table_NN$Week <- rownames(plot_table_NN)

plot_table_R <- rbind(pt_CC_ps_t_3_est, pt_CC_ps_t_2_est, 0, as.data.frame(Results_CC_ps$estimates))
colnames(plot_table_R) <- "estimate"
plot_table_R$high <- c(pt_CC_ps_t_3_high, pt_CC_ps_t_2_high, 0, Results_CC_ps$estimates + (Results_CC_ps$standard.error*1.96))
plot_table_R$low <- c(pt_CC_ps_t_3_low, pt_CC_ps_t_2_low, 0, Results_CC_ps$estimates - (Results_CC_ps$standard.error*1.96))
plot_table_R$Matching <- "Matching + Refinement" 
rownames(plot_table_R)[1:4] <- c("t-3","t-2","t-1","t")
plot_table_R$Week <- rownames(plot_table_R)

plot_table_CC <- rbind(plot_table_R, plot_table_NN)
plot_table_CC$Week <- factor(plot_table_CC$Week, levels = c("t-3", "t-2", "t-1", "t", "t+1", "t+2", "t+3", "t+4", "t+5"))

plot_table_CC$`Estimation Strategy` <- factor(plot_table_CC$Matching, levels = c("No Matching", "Matching + Refinement"))


# Plot
pd = position_dodge(.6)   
ggplot(plot_table_CC, aes(x = Week, y = estimate, color = `Estimation Strategy`)) + geom_hline(yintercept = 0, color = "gray40", linetype="dashed") +
  geom_point(shape = 16, size  = 6, position = pd) +
  geom_errorbar(aes(ymin  = low, ymax  = high), width = 0, size  = 1.45, position = pd) +
  theme_bw() + theme(axis.title = element_text(face = "bold")) + ylab("Rate of Change - New Weekly Cases Per Capita (IHST)") + scale_color_manual(values= c("gray60", "#588786")) + theme(legend.text=element_text(size=20), axis.text = element_text(size = 20), axis.title = element_text(size = 20), text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + ggtitle("Complete Closures") + scale_y_continuous(limits = c(-1.5, 1.5), breaks=c(-1.5,-1,-0.5,0,0.5,1,1.5))





# Analysis of Refinement strategies - change refinements method to compare

Complete_A <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "completenew_D", matching = TRUE, refinement.method = "CBPS.match", covs.formula = ~ human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older +  median_age + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + first.dose.pct + second.dose.pct + daily_tests_per_100k + I(lag(daily_tests_per_100k, 1:3))+ I(lag(Stringency_revised, 1:3)) + Stringency_revised + I(lag(partialnew_D, 1:3)) + partialnew_D + I(lag(Stringency_Treat, 1:3)) + Stringency_Treat, data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)


cvb_CC_A <- get_covariate_balance(Complete_A$att,DF_Full, covariates = c("gdp_per_capita_LOG", "new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "median_age", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy", "partialnew_D", "Stringency_Treat", "completenew_D"))

Var <- ncol(cvb_CC_A)
cvb_CC_A_totals <- data.frame(sum(ifelse(abs(cvb_CC_A[1,])>=0.20, abs(cvb_CC_A[1,])-0.20, 0))/Var, sum(ifelse(abs(cvb_CC_A[2,])>=0.20, abs(cvb_CC_A[2,])-0.20, 0))/Var, sum(ifelse(abs(cvb_CC_A[3,])>=0.20, abs(cvb_CC_A[3,])-0.20, 0))/Var, sum(ifelse(abs(cvb_CC_A[4,])>=0.20, abs(cvb_CC_A[4,])-0.20, 0))/Var)
cvb_CC_A_totals

# Imbalance scores by refinement strategy
  # None - 0.128
  # Mahal - 0.076
  # PS - 0.034
  # CBPS - 0.048


  # Partial Closure --------------------

Partial_NN <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", refinement.method = "none", treatment = "partialnew_D", matching = TRUE, data = DF_Full, match.missing = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

pt_PC_NN <- placeboTest(Partial_NN, data.in = DF_Full, plot = F)
pt_PC_NN_t_3_est <- pt_PC_NN$estimates[1]
pt_PC_NN_t_2_est <- pt_PC_NN$estimates[2]

pt_PC_NN_t_3_high <- quantile(pt_PC_NN$bootstrapped.estimates[ ,1],0.975)
pt_PC_NN_t_2_high <- quantile(pt_PC_NN$bootstrapped.estimates[ ,2],0.975)
pt_PC_NN_t_3_low <- quantile(pt_PC_NN$bootstrapped.estimates[ ,1],0.025)
pt_PC_NN_t_2_low <- quantile(pt_PC_NN$bootstrapped.estimates[ ,2],0.025)

Results_PC_NN <- PanelEstimate(sets = Partial_NN, data = DF_Full)



Partial_cbps <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "partialnew_D", refinement.method = "CBPS.match", covs.formula = ~ human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older +  median_age + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + first.dose.pct + second.dose.pct + daily_tests_per_100k + I(lag(daily_tests_per_100k, 1:3))+ I(lag(Stringency_revised, 1:3)) + Stringency_revised + I(lag(completenew_D, 1:3)) + completenew_D + I(lag(Stringency_Treat, 1:3)) + Stringency_Treat, data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

pt_PC_cbps <- placeboTest(Partial_cbps, data.in = DF_Full, plot = F)

pt_PC_cbps_t_3_est <- pt_PC_cbps$estimates[1]
pt_PC_cbps_t_2_est <- pt_PC_cbps$estimates[2]

pt_PC_cbps_t_3_high <- quantile(pt_PC_cbps$bootstrapped.estimates[ ,1],0.975)
pt_PC_cbps_t_2_high <- quantile(pt_PC_cbps$bootstrapped.estimates[ ,2],0.975)
pt_PC_cbps_t_3_low <- quantile(pt_PC_cbps$bootstrapped.estimates[ ,1],0.025)
pt_PC_cbps_t_2_low <- quantile(pt_PC_cbps$bootstrapped.estimates[ ,2],0.025)

Results_PC_cbps <- PanelEstimate(sets = Partial_cbps, data = DF_Full)
summary(Results_PC_cbps)






# Generate Figure 4

plot_table_NN <- rbind(pt_PC_NN_t_3_est, pt_PC_NN_t_2_est, 0, as.data.frame(Results_PC_NN$estimates))
colnames(plot_table_NN) <- "estimate"
plot_table_NN$high <- c(pt_PC_NN_t_3_high, pt_PC_NN_t_2_high, 0, Results_PC_NN$estimates + (Results_PC_NN$standard.error*1.96))
plot_table_NN$low <- c(pt_PC_NN_t_3_low, pt_PC_NN_t_2_low, 0, Results_PC_NN$estimates - (Results_PC_NN$standard.error*1.96))
plot_table_NN$Matching <- "No Matching" 
rownames(plot_table_NN)[1:4] <- c("t-3","t-2","t-1","t")
plot_table_NN$Week <- rownames(plot_table_NN)


plot_table_R <- rbind(pt_PC_cbps_t_3_est, pt_PC_cbps_t_2_est, 0, as.data.frame(Results_PC_cbps$estimates))
colnames(plot_table_R) <- "estimate"
plot_table_R$high <- c(pt_PC_cbps_t_3_high, pt_PC_cbps_t_2_high, 0, Results_PC_cbps$estimates + (Results_PC_cbps$standard.error*1.96))
plot_table_R$low <- c(pt_PC_cbps_t_3_low, pt_PC_cbps_t_2_low, 0, Results_PC_cbps$estimates - (Results_PC_cbps$standard.error*1.96))
plot_table_R$Matching <- "Matching + Refinement" 
rownames(plot_table_R)[1:4] <- c("t-3","t-2","t-1","t")
plot_table_R$Week <- rownames(plot_table_R)

plot_table_PC <- rbind(plot_table_R, plot_table_NN)
plot_table_PC$Week <- factor(plot_table_PC$Week, levels = c("t-3", "t-2", "t-1", "t", "t+1", "t+2", "t+3", "t+4", "t+5"))

plot_table_PC$`Estimation Strategy` <- factor(plot_table_PC$Matching, levels = c("No Matching", "Matching + Refinement"))



# Plot
pd = position_dodge(.6)   
ggplot(plot_table_PC, aes(x = Week, y = estimate, color = `Estimation Strategy`)) + geom_hline(yintercept = 0, color = "gray40", linetype="dashed") +
  geom_point(shape = 16, size  = 6, position = pd) +
  geom_errorbar(aes(ymin  = low, ymax  = high), width = 0, size  = 1.45, position = pd) +
  theme_bw() + theme(axis.title = element_text(face = "bold")) + ylab("Rate of Change - New Weekly Cases Per Capita (IHST)") + scale_color_manual(values= c("gray70", "#965453")) + theme(legend.text=element_text(size=20), axis.text = element_text(size = 20), axis.title = element_text(size = 20), text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + ggtitle("Partial Closures") + scale_y_continuous(limits = c(-1.5, 1.5), breaks=c(-1.5,-1,-0.5,0,0.5,1,1.5))







# Analysis of Refinement strategies

Partial_A <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "partialnew_D", matching = TRUE, refinement.method = "CBPS.match", covs.formula = ~ human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older +  median_age + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + first.dose.pct + second.dose.pct + daily_tests_per_100k + I(lag(daily_tests_per_100k, 1:3))+ I(lag(Stringency_revised, 1:3)) + Stringency_revised + I(lag(completenew_D, 1:3)) + completenew_D + I(lag(Stringency_Treat, 1:3)) + Stringency_Treat, data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)


cvb_PC_A <- get_covariate_balance(Partial_A$att,DF_Full, covariates = c("gdp_per_capita_LOG", "new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "median_age", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy", "partialnew_D", "Stringency_Treat", "completenew_D"))

Var <- ncol(cvb_PC_A)
cvb_PC_A_totals <- data.frame(sum(ifelse(abs(cvb_PC_A[1,])>=0.20, abs(cvb_PC_A[1,])-0.20, 0))/Var, sum(ifelse(abs(cvb_PC_A[2,])>=0.20, abs(cvb_PC_A[2,])-0.20, 0))/Var, sum(ifelse(abs(cvb_PC_A[3,])>=0.20, abs(cvb_PC_A[3,])-0.20, 0))/Var, sum(ifelse(abs(cvb_PC_A[4,])>=0.20, abs(cvb_PC_A[4,])-0.20, 0))/Var)
cvb_PC_A_totals


# Imbalance scores by refinement strategy
  # None - 0.246
  # Mahal - 0.068
  # PS - 0.009
  # CBPS - 0.006



  # Individual Closure Policies ----------------------------


# Analysis of Refinement strategies and Calculating T+2 estimates for Table 1. Replace Complete Closures sub-type in treatment.

Individual_A <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "essential_onlynew_D", refinement.method = "ps.match", matching =TRUE, covs.formula = ~ human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older + life_expectancy + gdp_per_capita_LOG + median_age + area_ifr_LOG  + first.dose.pct + second.dose.pct + Stringency_Treat + I(lag(Stringency_Treat, 1:3)) + daily_tests_per_100k + I(lag(daily_tests_per_100k, 1:3)) + Stringency_revised + I(lag(Stringency_revised, 1:3)) + specific_countrynew_D + I(lag(specific_countrynew_D, 1:3)) + work_excepnew_D + I(lag(work_excepnew_D, 1:3)) + citizen_excepnew_D + I(lag(citizen_excepnew_D, 1:3))  + essential_onlynew_D + I(lag(essential_onlynew_D, 1:3)) + partialnew_D + I(lag(partialnew_D, 1:3)) , data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = F)

#Results_ID_A <- PanelEstimate(sets = Individual_A, data = DF_Full)
#summary(Results_ID_A)

cvb_ID_A <- get_covariate_balance(Individual_A$att,DF_Full, covariates = c("gdp_per_capita_LOG", "new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "median_age", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy", "partialnew_D", "Stringency_Treat", "specific_countrynew_D", "work_excepnew_D", "citizen_excepnew_D", "essential_onlynew_D"))

Var <- ncol(cvb_ID_A)
cvb_ID_A_totals <- data.frame(sum(ifelse(abs(cvb_ID_A[1,])>=0.20, abs(cvb_ID_A[1,])-0.20, 0))/Var, sum(ifelse(abs(cvb_ID_A[2,])>=0.20, abs(cvb_ID_A[2,])-0.20, 0))/Var, sum(ifelse(abs(cvb_ID_A[3,])>=0.20, abs(cvb_ID_A[3,])-0.20, 0))/Var, sum(ifelse(abs(cvb_ID_A[4,])>=0.20, abs(cvb_ID_A[4,])-0.20, 0))/Var)
cvb_ID_A_totals


# Analysis of Refinement strategies

# specific_countrynew_D
  # none, T-1 = 0.556
  # Mahal, T-1 = 0.195
  # PS, T-1 = 0.0386
  # CBPS, T-1 = 0.0441

# Work Exception
  # none, T-1 = 0.293
  # Mahal, T-1 = 0.128
  # PS, T-1 = 0.021
  # CBPS, T-1 = 0.0098

# Citizen Exceptions
  # none, T-1 = 0.059
  # Mahal, T-1 = 0.058
  # PS, T-1 = 0.029
  # CBPS, T-1 = 0.049

# Essentials Only
  # none, T-1 = 0.109
  # Mahal, T-1 = 0.0387
  # PS, T-1 = 0.0431
  # CBPS, T-1 = 0.0282













Individual_B <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "visa_bannew_D", refinement.method = "CBPS.match", matching =TRUE, covs.formula = ~ human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older + daily_tests_per_100k + life_expectancy + gdp_per_capita_LOG + median_age + area_ifr_LOG  + Stringency_revised + first.dose.pct + second.dose.pct +  + Stringency_Treat + I(lag(Stringency_Treat, 1:3))  + I(lag(daily_tests_per_100k, 1:3)) + I(lag(Stringency_revised, 1:3)) + visa_bannew_D + I(lag(visa_bannew_D, 1:3)) + citizenship_bannew_D + I(lag(citizenship_bannew_D, 1:3)) + history_bannew_D + I(lag(history_bannew_D, 1:3))  + border_closurenew_D + I(lag(border_closurenew_D, 1:3)) + completenew_D + I(lag(completenew_D, 1:3)) , data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = F)

cvb_ID_B <- get_covariate_balance(Individual_B$att,DF_Full, covariates = c("gdp_per_capita_LOG", "new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "median_age", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy", "completenew_D", "Stringency_Treat", "visa_bannew_D", "citizenship_bannew_D", "history_bannew_D", "border_closurenew_D"))

Var <- ncol(cvb_ID_B)
cvb_ID_B_totals <- data.frame(sum(ifelse(abs(cvb_ID_B[1,])>=0.20, abs(cvb_ID_B[1,])-0.20, 0))/Var, sum(ifelse(abs(cvb_ID_B[2,])>=0.20, abs(cvb_ID_B[2,])-0.20, 0))/Var, sum(ifelse(abs(cvb_ID_B[3,])>=0.20, abs(cvb_ID_B[3,])-0.20, 0))/Var, sum(ifelse(abs(cvb_ID_B[4,])>=0.20, abs(cvb_ID_B[4,])-0.20, 0))/Var)
cvb_ID_B_totals


# Analysis of Refinement strategies

# Visa Ban
  # none, T-1 = 0.072
  # Mahal, T-1 = 0.051
  # PS, T-1 = 0.006
  # CBPS, T-1 = 0.001

# Citizenship Ban
  # none, T-1 = 0.368
  # Mahal, T-1 = 0.168
  # PS, T-1 = 0.029
  # CBPS, T-1 = 0.063

# Travel History 
  # None, T-1 = 0.244
  # Mahal, T-1 = 0.115
  # PS, T-1 = 0.0129
  # CBPS, T-1 = 0.0126

# Border Closure
  # None, T-1 = 0.127
  # Mahal, T-1 = 0.040
  # PS, T-1 = 0.002
  # CBPS, T-1 = 0.004


# Retrieve estimates for Table 1 - replace Treatment variable and refinement method accordingly

Individual_C <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "essential_onlynew_D", refinement.method = "CBPS.match", matching =TRUE, covs.formula = ~ human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older + daily_tests_per_100k + life_expectancy + gdp_per_capita_LOG + median_age + area_ifr_LOG   + first.dose.pct + second.dose.pct + Stringency_revised + I(lag(Stringency_revised, 1:3)) + Stringency_Treat + I(lag(Stringency_Treat, 1:3)) + I(lag(daily_tests_per_100k, 1:3)) + specific_countrynew_D + I(lag(specific_countrynew_D, 1:3)) + work_excepnew_D + I(lag(work_excepnew_D, 1:3)) + citizen_excepnew_D + I(lag(citizen_excepnew_D, 1:3))  + essential_onlynew_D + I(lag(essential_onlynew_D, 1:3)) + partialnew_D + I(lag(partialnew_D, 1:3)) , data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = F)

Results_Ind_C <- PanelEstimate(sets = Individual_C, data = DF_Full)
summary(Results_Ind_C)

Individual_C








Individual_D <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "border_closurenew_D", refinement.method = "ps.match", matching =TRUE, covs.formula = ~ human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older + daily_tests_per_100k + life_expectancy + gdp_per_capita_LOG + median_age + area_ifr_LOG  + Stringency_revised + first.dose.pct + second.dose.pct +  + Stringency_Treat + I(lag(Stringency_Treat, 1:3))  + I(lag(daily_tests_per_100k, 1:3)) + I(lag(Stringency_revised, 1:3)) + visa_bannew_D + I(lag(visa_bannew_D, 1:3)) + citizenship_bannew_D + I(lag(citizenship_bannew_D, 1:3)) + history_bannew_D + I(lag(history_bannew_D, 1:3))  + border_closurenew_D + I(lag(border_closurenew_D, 1:3)) + completenew_D + I(lag(completenew_D, 1:3)) , data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = F)

Individual_D

Results_Ind_D <- PanelEstimate(sets = Individual_D, data = DF_Full)
summary(Results_Ind_D)



  # Islands Only Analysis --------------

COBAP_Island_Countries_List <- read_excel("Islands_list.xlsx", col_names = F)

TEST_Islands <- subset(DF_Full, iso3c %in% COBAP_Island_Countries_List$...2)


# Complete Closures

Island_CC <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "completenew_D", refinement.method ="CBPS.match", matching = TRUE, covs.formula = ~ gdp_per_capita_LOG + area_ifr_LOG + median_age + first.dose.pct + second.dose.pct + life_expectancy + daily_tests_per_100k +  aged_70_older + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + Stringency_revised + I(lag(Stringency_Treat, 1:3)) + I(lag(daily_tests_per_100k, 1:3)) + partialnew_D + I(lag(partialnew_D, 1:3)), data = TEST_Islands, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = F)

cvb_IS_CC <- get_covariate_balance(Island_CC$att,DF_Full, covariates = c("gdp_per_capita_LOG", "new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "median_age", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy", "completenew_D", "Stringency_Treat", "partialnew_D"))

Var <- ncol(cvb_IS_CC)
cvb_IS_CC_totals <- data.frame(sum(ifelse(abs(cvb_IS_CC[1,])>=0.20, abs(cvb_IS_CC[1,])-0.20, 0))/Var, sum(ifelse(abs(cvb_IS_CC[2,])>=0.20, abs(cvb_IS_CC[2,])-0.20, 0))/Var, sum(ifelse(abs(cvb_IS_CC[3,])>=0.20, abs(cvb_IS_CC[3,])-0.20, 0))/Var, sum(ifelse(abs(cvb_IS_CC[4,])>=0.20, abs(cvb_IS_CC[4,])-0.20, 0))/Var)
cvb_IS_CC_totals


# None - 0.128
# Mahal - 0.107
# PS - 0.116
# CBPS - 0.112





Island_CC_NN <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "completenew_D", refinement.method ="none", matching = FALSE, covs.formula = ~ gdp_per_capita_LOG + area_ifr_LOG + median_age + first.dose.pct + second.dose.pct + life_expectancy + daily_tests_per_100k +  aged_70_older + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG +  Stringency_revised + I(lag(Stringency_Treat, 1:3)) + I(lag(daily_tests_per_100k, 1:3)) + partialnew_D + I(lag(partialnew_D, 1:3)), data = TEST_Islands, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = F)

Results_IS_CC_NN <- PanelEstimate(sets = Island_CC_NN, data = TEST_Islands)

pt_IS_CC_NN <- placeboTest(Island_CC_NN, data.in = TEST_Islands, plot = F)
pt_IS_CC_NN_t_3_est <- pt_IS_CC_NN$estimates[1]
pt_IS_CC_NN_t_2_est <- pt_IS_CC_NN$estimates[2]

pt_IS_CC_NN_t_3_high <- quantile(pt_IS_CC_NN$bootstrapped.estimates[ ,1],0.975)
pt_IS_CC_NN_t_2_high <- quantile(pt_IS_CC_NN$bootstrapped.estimates[ ,2],0.975)
pt_IS_CC_NN_t_3_low <- quantile(pt_IS_CC_NN$bootstrapped.estimates[ ,1],0.025)
pt_IS_CC_NN_t_2_low <- quantile(pt_IS_CC_NN$bootstrapped.estimates[ ,2],0.025)



Island_CC_m <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "completenew_D", refinement.method ="none", matching = TRUE, covs.formula = ~ gdp_per_capita_LOG + area_ifr_LOG + median_age + first.dose.pct + second.dose.pct + life_expectancy + daily_tests_per_100k +  aged_70_older + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG +  Stringency_revised + I(lag(Stringency_Treat, 1:3)) + I(lag(daily_tests_per_100k, 1:3)) + partialnew_D + I(lag(partialnew_D, 1:3)), data = TEST_Islands, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = F)

Results_IS_CC_m <- PanelEstimate(sets = Island_CC_m, data = TEST_Islands)
summary(Results_IS_CC_m)

pt_IS_CC_m <- placeboTest(Island_CC_m, data.in = TEST_Islands, plot = F)
pt_IS_CC_m_t_3_est <- pt_IS_CC_m$estimates[1]
pt_IS_CC_m_t_2_est <- pt_IS_CC_m$estimates[2]

pt_IS_CC_m_t_3_high <- quantile(pt_IS_CC_m$bootstrapped.estimates[ ,1],0.975)
pt_IS_CC_m_t_2_high <- quantile(pt_IS_CC_m$bootstrapped.estimates[ ,2],0.975)
pt_IS_CC_m_t_3_low <- quantile(pt_IS_CC_m$bootstrapped.estimates[ ,1],0.025)
pt_IS_CC_m_t_2_low <- quantile(pt_IS_CC_m$bootstrapped.estimates[ ,2],0.025)





# Generate Figure 5

plot_table_NN <- rbind(pt_IS_CC_NN_t_3_est, pt_IS_CC_NN_t_2_est, 0, as.data.frame(Results_Is_CC_NN$estimates))
colnames(plot_table_NN) <- "estimate"
plot_table_NN$high <- c(pt_IS_CC_NN_t_3_high, pt_IS_CC_NN_t_2_high, 0, Results_Is_CC_NN$estimates + (Results_Is_CC_NN$standard.error*1.96))
plot_table_NN$low <- c(pt_IS_CC_NN_t_3_low, pt_IS_CC_NN_t_2_low, 0, Results_Is_CC_NN$estimates - (Results_Is_CC_NN$standard.error*1.96))
plot_table_NN$Matching <- "No Matching" 
rownames(plot_table_NN)[1:4] <- c("t-3","t-2","t-1","t")
plot_table_NN$Week <- rownames(plot_table_NN)


plot_table_m <- rbind(pt_IS_CC_m_t_3_est, pt_IS_CC_m_t_2_est, 0, as.data.frame(Results_Is_CC_m$estimates))
colnames(plot_table_m) <- "estimate"
plot_table_m$high <- c(pt_IS_CC_m_t_3_high, pt_IS_CC_m_t_2_high, 0, Results_Is_CC_m$estimates + (Results_Is_CC_m$standard.error*1.96))
plot_table_m$low <- c(pt_IS_CC_m_t_3_low, pt_IS_CC_m_t_2_low, 0, Results_Is_CC_m$estimates - (Results_Is_CC_m$standard.error*1.96))
plot_table_m$Matching <- "Matching" 
rownames(plot_table_m)[1:4] <- c("t-3","t-2","t-1","t")
plot_table_m$Week <- rownames(plot_table_m)



plot_table_CC <- rbind(plot_table_NN, plot_table_m)
plot_table_CC$Week <- factor(plot_table_CC$Week, levels = c("t-3", "t-2", "t-1", "t", "t+1", "t+2", "t+3", "t+4", "t+5"))

plot_table_CC$`Estimation Strategy` <- factor(plot_table_CC$Matching, levels = c("No Matching", "Matching"))



# Plot
pd = position_dodge(.6)   
ggplot(plot_table_CC, aes(x = Week, y = estimate, color = `Estimation Strategy`)) + geom_hline(yintercept = 0, color = "gray40", linetype="dashed") +
  geom_point(shape = 16, size  = 6, position = pd) +
  geom_errorbar(aes(ymin  = low, ymax  = high), width = 0, size  = 1.45, position = pd) +
  theme_bw() + theme(axis.title = element_text(face = "bold")) + ylab("Rate of Change - New Weekly Cases Per Capita (IHST)") + scale_color_manual(values= c("gray70", "#588786")) + theme(legend.text=element_text(size=20), axis.text = element_text(size = 20), axis.title = element_text(size = 20), text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + ggtitle("Islands subset - Complete Closures") + scale_y_continuous(limits = c(-2.0, 2.5), breaks=c(-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5))









# Partial Closures

Island_PC <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "partialnew_D", refinement.method ="ps.match", matching = TRUE, covs.formula = ~ gdp_per_capita_LOG + area_ifr_LOG + median_age + first.dose.pct + second.dose.pct + life_expectancy +  daily_tests_per_100k +  aged_70_older + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG +  Stringency_revised + Stringency_Treat + I(lag(Stringency_Treat, 1:3)) + I(lag(daily_tests_per_100k, 1:3)) + completenew_D + I(lag(completenew_D, 1:3)), data = TEST_Islands, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)


cvb_IS_PC <- get_covariate_balance(Island_PC$att,DF_Full, covariates = c("gdp_per_capita_LOG", "new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "median_age", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy", "completenew_D", "Stringency_Treat", "partialnew_D"))

Var <- ncol(cvb_IS_PC)
cvb_IS_PC_totals <- data.frame(sum(ifelse(abs(cvb_IS_PC[1,])>=0.20, abs(cvb_IS_PC[1,])-0.20, 0))/Var, sum(ifelse(abs(cvb_IS_PC[2,])>=0.20, abs(cvb_IS_PC[2,])-0.20, 0))/Var, sum(ifelse(abs(cvb_IS_PC[3,])>=0.20, abs(cvb_IS_PC[3,])-0.20, 0))/Var, sum(ifelse(abs(cvb_IS_PC[4,])>=0.20, abs(cvb_IS_PC[4,])-0.20, 0))/Var)
cvb_IS_PC_totals


# None - 0.364
# Mahal - 0.334
# PS - 0.289
# CBPS - 0.296







Island_PC_NN <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "partialnew_D", refinement.method ="none", matching = FALSE, covs.formula = ~ gdp_per_capita_LOG + area_ifr_LOG + median_age + first.dose.pct + second.dose.pct + life_expectancy +  daily_tests_per_100k +  aged_70_older + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG +  Stringency_revised + Stringency_Treat + I(lag(Stringency_Treat, 1:3)) + I(lag(daily_tests_per_100k, 1:3)) + completenew_D + I(lag(completenew_D, 1:3)), data = TEST_Islands, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

Results_Is_PC_NN <- PanelEstimate(sets = Island_PC_NN, data = TEST_Islands)

pt_Is_PC_NN <- placeboTest(Island_PC_NN, data.in = TEST_Islands, plot = F)
pt_Is_PC_NN_t_3_est <- pt_Is_PC_NN$estimates[1]
pt_Is_PC_NN_t_2_est <- pt_Is_PC_NN$estimates[2]

pt_Is_PC_NN_t_3_high <- quantile(pt_Is_PC_NN$bootstrapped.estimates[ ,1],0.975)
pt_Is_PC_NN_t_2_high <- quantile(pt_Is_PC_NN$bootstrapped.estimates[ ,2],0.975)
pt_Is_PC_NN_t_3_low <- quantile(pt_Is_PC_NN$bootstrapped.estimates[ ,1],0.025)
pt_Is_PC_NN_t_2_low <- quantile(pt_Is_PC_NN$bootstrapped.estimates[ ,2],0.025)



Island_PC_ps <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "partialnew_D", refinement.method ="ps.match", matching = TRUE, covs.formula = ~ gdp_per_capita_LOG + area_ifr_LOG + median_age + first.dose.pct + second.dose.pct + life_expectancy +  daily_tests_per_100k +  aged_70_older + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG +  Stringency_revised + I(lag(Stringency_Treat, 1:3)) + I(lag(daily_tests_per_100k, 1:3)) + completenew_D + I(lag(completenew_D, 1:3)), data = TEST_Islands, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

Results_Is_PC_ps <- PanelEstimate(sets = Island_PC_ps, data = TEST_Islands)

pt_Is_PC_ps <- placeboTest(Island_PC_ps, data.in = TEST_Islands, plot = F)
pt_Is_PC_ps_t_3_est <- pt_Is_PC_ps$estimates[1]
pt_Is_PC_ps_t_2_est <- pt_Is_PC_ps$estimates[2]

pt_Is_PC_ps_t_3_high <- quantile(pt_Is_PC_ps$bootstrapped.estimates[ ,1],0.975)
pt_Is_PC_ps_t_2_high <- quantile(pt_Is_PC_ps$bootstrapped.estimates[ ,2],0.975)
pt_Is_PC_ps_t_3_low <- quantile(pt_Is_PC_ps$bootstrapped.estimates[ ,1],0.025)
pt_Is_PC_ps_t_2_low <- quantile(pt_Is_PC_ps$bootstrapped.estimates[ ,2],0.025)




# Generate Figure

plot_table_NN <- rbind(pt_Is_PC_NN_t_3_est, pt_Is_PC_NN_t_2_est, 0, as.data.frame(Results_Is_PC_NN$estimates))
colnames(plot_table_NN) <- "estimate"
plot_table_NN$high <- c(pt_Is_PC_NN_t_3_high, pt_Is_PC_NN_t_2_high, 0, Results_Is_PC_NN$estimates + (Results_Is_PC_NN$standard.error*1.96))
plot_table_NN$low <- c(pt_Is_PC_NN_t_3_low, pt_Is_PC_NN_t_2_low, 0, Results_Is_PC_NN$estimates - (Results_Is_PC_NN$standard.error*1.96))
plot_table_NN$Matching <- "No Matching" 
rownames(plot_table_NN)[1:4] <- c("t-3","t-2","t-1","t")
plot_table_NN$Week <- rownames(plot_table_NN)

plot_table_ps <- rbind(pt_Is_PC_ps_t_3_est, pt_Is_PC_ps_t_2_est, 0, as.data.frame(Results_Is_PC_ps$estimates))
colnames(plot_table_ps) <- "estimate"
plot_table_ps$high <- c(pt_Is_PC_ps_t_3_high, pt_Is_PC_ps_t_2_high, 0, Results_Is_PC_ps$estimates + (Results_Is_PC_ps$standard.error*1.96))
plot_table_ps$low <- c(pt_Is_PC_ps_t_3_low, pt_Is_PC_ps_t_2_low, 0, Results_Is_PC_ps$estimates - (Results_Is_PC_ps$standard.error*1.96))
plot_table_ps$Matching <- "Matching + Refinement" 
rownames(plot_table_ps)[1:4] <- c("t-3","t-2","t-1","t")
plot_table_ps$Week <- rownames(plot_table_ps)

plot_table_PC <- rbind(plot_table_NN, plot_table_ps)
plot_table_PC$Week <- factor(plot_table_PC$Week, levels = c("t-3", "t-2", "t-1", "t", "t+1", "t+2", "t+3", "t+4", "t+5"))

plot_table_PC$`Estimation Strategy` <- factor(plot_table_PC$Matching, levels = c("No Matching", "Matching + Refinement"))



# Plot
pd = position_dodge(.6)   
ggplot(plot_table_PC, aes(x = Week, y = estimate, color = `Estimation Strategy`)) + geom_hline(yintercept = 0, color = "gray40", linetype="dashed") +
  geom_point(shape = 16, size  = 6, position = pd) +
  geom_errorbar(aes(ymin  = low, ymax  = high), width = 0, size  = 1.45, position = pd) +
  theme_bw() + theme(axis.title = element_text(face = "bold")) + ylab("Rate of Change - New Weekly Cases Per Capita (IHST)") + scale_color_manual(values= c("gray70", "#965453")) + theme(legend.text=element_text(size=20), axis.text = element_text(size = 20), axis.title = element_text(size = 20), text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + ggtitle("Islands subset - Partial Closures") + scale_y_continuous(limits = c(-2.0, 2.5), breaks=c(-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5))










# Lockdowns

Island_DL <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "Stringency_Treat", refinement.method ="CBPS.match", matching = TRUE, covs.formula = ~ gdp_per_capita_LOG + area_ifr_LOG + median_age + first.dose.pct + second.dose.pct + life_expectancy +  daily_tests_per_100k +  aged_70_older + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + Stringency_revised +  partialnew_D + I(lag(partialnew_D, 1:3)) + I(lag(daily_tests_per_100k, 1:3)) + completenew_D + I(lag(completenew_D, 1:3)), data = TEST_Islands, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)


cvb_Is_DL <- get_covariate_balance(Island_DL$att,TEST_Islands, covariates = c("gdp_per_capita_LOG", "new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy", "median_age", "completenew_D", "partialnew_D"))

Var <- ncol(cvb_Is_DL)
cvb_IS_PC_totals <- data.frame(sum(ifelse(abs(cvb_Is_DL[1,])>=0.20, abs(cvb_Is_DL[1,])-0.20, 0))/Var, sum(ifelse(abs(cvb_Is_DL[2,])>=0.20, abs(cvb_Is_DL[2,])-0.20, 0))/Var, sum(ifelse(abs(cvb_Is_DL[3,])>=0.20, abs(cvb_Is_DL[3,])-0.20, 0))/Var, sum(ifelse(abs(cvb_Is_DL[4,])>=0.20, abs(cvb_Is_DL[4,])-0.20, 0))/Var)
cvb_IS_PC_totals

# Islands DL
  # None = 0.188
  # Mahal = 0.110
  # PS = 0.143
  # CBPS = 0.124



Island_DL_NN <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "Stringency_Treat", refinement.method ="none", matching = FALSE, covs.formula = ~ gdp_per_capita_LOG + area_ifr_LOG + median_age + first.dose.pct + second.dose.pct + life_expectancy +  daily_tests_per_100k +  aged_70_older + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG +  Stringency_revised + I(lag(partialnew_D, 1:3)) + I(lag(daily_tests_per_100k, 1:3)) + completenew_D + I(lag(completenew_D, 1:3)), data = TEST_Islands, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

Results_Is_DL_NN <- PanelEstimate(sets = Island_DL_NN, data = TEST_Islands)

pt_Is_DL_NN <- placeboTest(Island_DL_NN, data.in = TEST_Islands, plot = F)
pt_Is_DL_NN_t_3_est <- pt_Is_DL_NN$estimates[1]
pt_Is_DL_NN_t_2_est <- pt_Is_DL_NN$estimates[2]

pt_Is_DL_NN_t_3_high <- quantile(pt_Is_DL_NN$bootstrapped.estimates[ ,1],0.975)
pt_Is_DL_NN_t_2_high <- quantile(pt_Is_DL_NN$bootstrapped.estimates[ ,2],0.975)
pt_Is_DL_NN_t_3_low <- quantile(pt_Is_DL_NN$bootstrapped.estimates[ ,1],0.025)
pt_Is_DL_NN_t_2_low <- quantile(pt_Is_DL_NN$bootstrapped.estimates[ ,2],0.025)



Island_DL_ps <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "Stringency_Treat", refinement.method ="ps.match", matching = TRUE, covs.formula = ~ gdp_per_capita_LOG + area_ifr_LOG + median_age + first.dose.pct + second.dose.pct + life_expectancy +  daily_tests_per_100k +  aged_70_older + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + partialnew_D + Stringency_revised + I(lag(partialnew_D, 1:3)) + I(lag(daily_tests_per_100k, 1:3)) + completenew_D + I(lag(completenew_D, 1:3)), data = TEST_Islands, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

Results_Is_DL_ps <- PanelEstimate(sets = Island_DL_ps, data = TEST_Islands)

pt_Is_DL_ps <- placeboTest(Island_DL_ps, data.in = TEST_Islands, plot = F)
pt_Is_DL_ps_t_3_est <- pt_Is_DL_ps$estimates[1]
pt_Is_DL_ps_t_2_est <- pt_Is_DL_ps$estimates[2]

pt_Is_DL_ps_t_3_high <- quantile(pt_Is_DL_ps$bootstrapped.estimates[ ,1],0.975)
pt_Is_DL_ps_t_2_high <- quantile(pt_Is_DL_ps$bootstrapped.estimates[ ,2],0.975)
pt_Is_DL_ps_t_3_low <- quantile(pt_Is_DL_ps$bootstrapped.estimates[ ,1],0.025)
pt_Is_DL_ps_t_2_low <- quantile(pt_Is_DL_ps$bootstrapped.estimates[ ,2],0.025)






# Generate Figure

plot_table_NN <- rbind(pt_Is_DL_NN_t_3_est, pt_Is_DL_NN_t_2_est, 0, as.data.frame(Results_Is_DL_NN$estimates))
colnames(plot_table_NN) <- "estimate"
plot_table_NN$high <- c(pt_Is_DL_NN_t_3_high, pt_Is_DL_NN_t_2_high, 0, Results_Is_DL_NN$estimates + (Results_Is_DL_NN$standard.error*1.96))
plot_table_NN$low <- c(pt_Is_DL_NN_t_3_low, pt_Is_DL_NN_t_2_low, 0, Results_Is_DL_NN$estimates - (Results_Is_DL_NN$standard.error*1.96))
plot_table_NN$Matching <- "No Matching" 
rownames(plot_table_NN)[1:4] <- c("t-3","t-2","t-1","t")
plot_table_NN$Week <- rownames(plot_table_NN)

plot_table_ps <- rbind(pt_Is_DL_ps_t_3_est, pt_Is_DL_ps_t_2_est, 0, as.data.frame(Results_Is_DL_ps$estimates))
colnames(plot_table_ps) <- "estimate"
plot_table_ps$high <- c(pt_Is_DL_ps_t_3_high, pt_Is_DL_ps_t_2_high, 0, Results_Is_DL_ps$estimates + (Results_Is_DL_ps$standard.error*1.96))
plot_table_ps$low <- c(pt_Is_DL_ps_t_3_low, pt_Is_DL_ps_t_2_low, 0, Results_Is_DL_ps$estimates - (Results_Is_DL_ps$standard.error*1.96))
plot_table_ps$Matching <- "Matching + Refinement" 
rownames(plot_table_ps)[1:4] <- c("t-3","t-2","t-1","t")
plot_table_ps$Week <- rownames(plot_table_ps)


plot_table_DL <- rbind(plot_table_NN, plot_table_ps)
plot_table_DL$Week <- factor(plot_table_DL$Week, levels = c("t-3", "t-2", "t-1", "t", "t+1", "t+2", "t+3", "t+4", "t+5"))

plot_table_DL$`Estimation Strategy` <- factor(plot_table_DL$Matching, levels = c("No Matching", "Matching + Refinement"))



# Plot
pd = position_dodge(.6)   
ggplot(plot_table_DL, aes(x = Week, y = estimate, color = `Estimation Strategy`)) + geom_hline(yintercept = 0, color = "gray40", linetype="dashed") +
  geom_point(shape = 16, size  = 6, position = pd) +
  geom_errorbar(aes(ymin  = low, ymax  = high), width = 0, size  = 1.45, position = pd) +
  theme_bw() + theme(axis.title = element_text(face = "bold")) + ylab("Rate of Change - New Weekly Cases Per Capita (IHST)") + scale_color_manual(values= c("gray70", "#689950")) + theme(legend.text=element_text(size=20), axis.text = element_text(size = 20), axis.title = element_text(size = 20), text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + ggtitle("Islands subset - Domestic Lockdowns") + scale_y_continuous(limits = c(-2.0, 2.5), breaks=c(-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5))



  # Covariate Balancing Plots ----------

# Lockdowns

Lockdown_N <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", matching=FALSE, treatment = "Stringency_Treat", refinement.method = "none",  covs.formula = ~ human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older +  median_age + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + first.dose.pct + second.dose.pct + daily_tests_per_100k + I(lag(daily_tests_per_100k, 1:3)) + Stringency_revised + I(lag(completenew_D, 1:3)) + completenew_D + I(lag(partialnew_D, 1:3)) + partialnew_D, data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

cvb_LD_1A <- as.data.frame(get_covariate_balance(Lockdown_N$att,DF_Full, covariates = c("human_development_index", "v2x_libdem", "v2x_freexp_altinf", "pop.2020_LOG", "aged_70_older", "median_age", "life_expectancy",  "gdp_per_capita_LOG",  "area_ifr_LOG")))

cvb_LD_1B <- as.data.frame(get_covariate_balance(Lockdown_N$att,DF_Full, covariates = c("new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "partialnew_D", "completenew_D")))



Lockdown_m <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "Stringency_Treat", refinement.method = "mahalanobis", covs.formula = ~ human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older +  median_age + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + first.dose.pct + second.dose.pct + daily_tests_per_100k + I(lag(daily_tests_per_100k, 1:3)) + Stringency_revised + I(lag(completenew_D, 1:3)) + completenew_D + I(lag(partialnew_D, 1:3)) + partialnew_D, data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

cvb_LD_2A <- as.data.frame(get_covariate_balance(Lockdown_m$att,DF_Full, covariates = c("median_age", "gdp_per_capita_LOG", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy")))

cvb_LD_2B <- as.data.frame(get_covariate_balance(Lockdown_m$att,DF_Full, covariates = c("new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "partialnew_D", "completenew_D")))




Lockdown_ps <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "Stringency_Treat", refinement.method = "ps.match", covs.formula = ~ human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older +  median_age + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + first.dose.pct + second.dose.pct + daily_tests_per_100k + I(lag(daily_tests_per_100k, 1:3)) + Stringency_revised + I(lag(completenew_D, 1:3)) + completenew_D + I(lag(partialnew_D, 1:3)) + partialnew_D, data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

cvb_LD_3A <- as.data.frame(get_covariate_balance(Lockdown_ps$att,DF_Full, covariates = c("median_age", "gdp_per_capita_LOG", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy")))

cvb_LD_3B <- as.data.frame(get_covariate_balance(Lockdown_ps$att,DF_Full, covariates = c("new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "partialnew_D", "completenew_D")))



Lockdown_cbps <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "Stringency_Treat", refinement.method = "CBPS.match", covs.formula = ~ human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older +  median_age + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + first.dose.pct + second.dose.pct + daily_tests_per_100k + I(lag(daily_tests_per_100k, 1:3)) + Stringency_revised + I(lag(completenew_D, 1:3)) + completenew_D + I(lag(partialnew_D, 1:3)) + partialnew_D, data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

cvb_LD_4A <- as.data.frame(get_covariate_balance(Lockdown_cbps$att,DF_Full, covariates = c("median_age", "gdp_per_capita_LOG", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy")))

cvb_LD_4B <- as.data.frame(get_covariate_balance(Lockdown_cbps$att,DF_Full, covariates = c("new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "partialnew_D", "completenew_D")))


cvb_LD_1A <-apply(as.data.frame(apply(cvb_LD_1A, 1, abs)),2, mean)
cvb_LD_1B <-apply(as.data.frame(apply(cvb_LD_1B, 1, abs)),2, mean)
cvb_LD_2A <-apply(as.data.frame(apply(cvb_LD_2A, 1, abs)),2, mean)
cvb_LD_2B <-apply(as.data.frame(apply(cvb_LD_2B, 1, abs)),2, mean)
cvb_LD_3A <-apply(as.data.frame(apply(cvb_LD_3A, 1, abs)),2, mean)
cvb_LD_3B <-apply(as.data.frame(apply(cvb_LD_3B, 1, abs)),2, mean)
cvb_LD_4A <-apply(as.data.frame(apply(cvb_LD_4A, 1, abs)),2, mean)
cvb_LD_4B <-apply(as.data.frame(apply(cvb_LD_4B, 1, abs)),2, mean)


cvb_plot <- as.data.frame(rbind(cvb_LD_1A, cvb_LD_1B, cvb_LD_2A, cvb_LD_2B, cvb_LD_3A, cvb_LD_3B, cvb_LD_4A, cvb_LD_4B))

cvb_plot$Refinement <- factor(c("None","None","Mahalanobis","Mahalanobis","PS","PS","CBPS","CBPS"), levels= c("None", "Mahalanobis", "PS", "CBPS"))

cvb_plot$`Variable Type` <- as.factor(c("Fixed","Time-Variant","Fixed","Time-Variant","Fixed","Time-Variant","Fixed","Time-Variant"))

cvb_plot <- melt(cvb_plot, id.vars=c("Refinement", "Variable Type"))
colnames(cvb_plot)[3] <- "Time Period"
colnames(cvb_plot)[4] <- "Average of Standardized Mean Differences"

ggplot(cvb_plot, aes(`Time Period`, `Average of Standardized Mean Differences`, group = interaction(Refinement, `Variable Type`),  color = Refinement, linetype = `Variable Type`)) + geom_hline(yintercept=0.20, color = "black", size=0.9)+ geom_line(size=1.8) + geom_point(size=6.5, shape=18) + theme_bw(base_size = 22) + scale_y_continuous(breaks=seq(0.0, 0.6, 0.1), limits=c(0, 0.6))









# Complete

Complete_N <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "completenew_D", refinement.method = "none", covs.formula = ~ human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older + median_age + daily_tests_per_100k + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + partialnew_D + Stringency_Treat + Stringency_revised + first.dose.pct + second.dose.pct + I(lag(Stringency_revised, 1:3)) + I(lag(partialnew_D, 1:3)) + I(lag(Stringency_Treat, 1:3))  + I(lag(daily_tests_per_100k, 1:3)), data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

cvb_CC_1A <- as.data.frame(get_covariate_balance(Complete_N$att,DF_Full, covariates = c("median_age", "gdp_per_capita_LOG", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy")))

cvb_CC_1B <- as.data.frame(get_covariate_balance(Complete_N$att,DF_Full, covariates = c("new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "Stringency_Treat", "partialnew_D")))



Complete_m <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "completenew_D", refinement.method = "mahalanobis", covs.formula = ~ median_age + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older + daily_tests_per_100k + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + partialnew_D + Stringency_Treat + Stringency_revised + first.dose.pct + second.dose.pct + I(lag(Stringency_revised, 1:3)) + I(lag(partialnew_D, 1:3)) + I(lag(Stringency_Treat, 1:3))  + I(lag(daily_tests_per_100k, 1:3)), data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

cvb_CC_2A <- as.data.frame(get_covariate_balance(Complete_m$att,DF_Full, covariates = c("median_age", "gdp_per_capita_LOG", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy")))

cvb_CC_2B <- as.data.frame(get_covariate_balance(Complete_m$att,DF_Full, covariates = c("new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "Stringency_Treat", "partialnew_D")))




Complete_ps <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "completenew_D", refinement.method = "ps.match", covs.formula = ~ median_age + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older + daily_tests_per_100k + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + partialnew_D + Stringency_Treat + Stringency_revised + first.dose.pct + second.dose.pct + I(lag(Stringency_revised, 1:3)) + I(lag(partialnew_D, 1:3)) + I(lag(Stringency_Treat, 1:3))  + I(lag(daily_tests_per_100k, 1:3)), data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

cvb_CC_3A <- as.data.frame(get_covariate_balance(Complete_ps$att,DF_Full, covariates = c("median_age", "gdp_per_capita_LOG", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy")))

cvb_CC_3B <- as.data.frame(get_covariate_balance(Complete_ps$att,DF_Full, covariates = c("new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "Stringency_Treat", "partialnew_D")))



Complete_cbps <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "completenew_D", refinement.method = "CBPS.match", covs.formula = ~ median_age + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older + daily_tests_per_100k + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + partialnew_D + Stringency_Treat + Stringency_revised + first.dose.pct + second.dose.pct + I(lag(Stringency_revised, 1:3)) + I(lag(partialnew_D, 1:3)) + I(lag(Stringency_Treat, 1:3))  + I(lag(daily_tests_per_100k, 1:3)), data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

cvb_CC_4A <- as.data.frame(get_covariate_balance(Complete_cbps$att,DF_Full, covariates = c("median_age", "gdp_per_capita_LOG", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy")))

cvb_CC_4B <- as.data.frame(get_covariate_balance(Complete_cbps$att,DF_Full, covariates = c("new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "Stringency_Treat", "partialnew_D")))


cvb_CC_1A <-apply(as.data.frame(apply(cvb_CC_1A, 1, abs)),2, mean)
cvb_CC_1B <-apply(as.data.frame(apply(cvb_CC_1B, 1, abs)),2, mean)
cvb_CC_2A <-apply(as.data.frame(apply(cvb_CC_2A, 1, abs)),2, mean)
cvb_CC_2B <-apply(as.data.frame(apply(cvb_CC_2B, 1, abs)),2, mean)
cvb_CC_3A <-apply(as.data.frame(apply(cvb_CC_3A, 1, abs)),2, mean)
cvb_CC_3B <-apply(as.data.frame(apply(cvb_CC_3B, 1, abs)),2, mean)
cvb_CC_4A <-apply(as.data.frame(apply(cvb_CC_4A, 1, abs)),2, mean)
cvb_CC_4B <-apply(as.data.frame(apply(cvb_CC_4B, 1, abs)),2, mean)


cvb_plot <- as.data.frame(rbind(cvb_CC_1A, cvb_CC_1B, cvb_CC_2A, cvb_CC_2B, cvb_CC_3A, cvb_CC_3B, cvb_CC_4A, cvb_CC_4B))

cvb_plot$Refinement <- factor(c("None","None","Mahalanobis","Mahalanobis","PS","PS","CBPS","CBPS"), levels= c("None", "Mahalanobis", "PS", "CBPS"))

cvb_plot$`Variable Type` <- as.factor(c("Fixed","Time-Variant","Fixed","Time-Variant","Fixed","Time-Variant","Fixed","Time-Variant"))


cvb_plot <- melt(cvb_plot, id.vars=c("Refinement", "Variable Type"))
colnames(cvb_plot)[3] <- "Time Period"
colnames(cvb_plot)[4] <- "Average of Standardized Mean Differences"

ggplot(cvb_plot, aes(`Time Period`, `Average of Standardized Mean Differences`, group = interaction(Refinement, `Variable Type`),  color = Refinement, linetype = `Variable Type`)) + geom_hline(yintercept=0.20, color = "black", size=0.9)+ geom_line(size=1.8) + geom_point(size=6.5, shape=18) + theme_bw(base_size = 22) + scale_y_continuous(breaks=seq(0.0, 0.6, 0.1), limits=c(0, 0.6))









# Partial

Partial_N <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "partialnew_D", refinement.method = "none", covs.formula = ~ median_age + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older + daily_tests_per_100k + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + completenew_D + Stringency_Treat + Stringency_revised + first.dose.pct + second.dose.pct + I(lag(Stringency_revised, 1:3)) + I(lag(completenew_D, 1:3)) + I(lag(Stringency_Treat, 1:3))  + I(lag(daily_tests_per_100k, 1:3)), data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

cvb_PC_1A <- as.data.frame(get_covariate_balance(Partial_N$att,DF_Full, covariates = c("median_age", "gdp_per_capita_LOG", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy")))

cvb_PC_1B <- as.data.frame(get_covariate_balance(Partial_N$att,DF_Full, covariates = c("new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "Stringency_Treat", "completenew_D")))



Partial_m <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "partialnew_D", refinement.method = "mahalanobis", covs.formula = ~ median_age + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older + daily_tests_per_100k + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + completenew_D + Stringency_Treat + Stringency_revised + first.dose.pct + second.dose.pct + I(lag(Stringency_revised, 1:3)) + I(lag(completenew_D, 1:3)) + I(lag(Stringency_Treat, 1:3))  + I(lag(daily_tests_per_100k, 1:3)), data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

cvb_PC_2A <- as.data.frame(get_covariate_balance(Partial_m$att,DF_Full, covariates = c("median_age", "gdp_per_capita_LOG", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy")))

cvb_PC_2B <- as.data.frame(get_covariate_balance(Partial_m$att,DF_Full, covariates = c("median_age", "new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "Stringency_Treat", "completenew_D")))




Partial_ps <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "partialnew_D", refinement.method = "ps.match", covs.formula = ~ median_age + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older + daily_tests_per_100k + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + completenew_D + Stringency_Treat + Stringency_revised + first.dose.pct + second.dose.pct + I(lag(Stringency_revised, 1:3)) + I(lag(completenew_D, 1:3)) + I(lag(Stringency_Treat, 1:3))  + I(lag(daily_tests_per_100k, 1:3)), data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

cvb_PC_3A <- as.data.frame(get_covariate_balance(Partial_ps$att,DF_Full, covariates = c("median_age", "gdp_per_capita_LOG", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy")))

cvb_PC_3B <- as.data.frame(get_covariate_balance(Partial_ps$att,DF_Full, covariates = c("median_age", "new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "Stringency_Treat", "completenew_D")))



Partial_cbps <- PanelMatch(lag = 3, time.id = "week_count", unit.id = "countryN", treatment = "partialnew_D", refinement.method = "CBPS.match", covs.formula = ~ median_age + human_development_index + v2x_libdem + v2x_freexp_altinf + pop.2020_LOG + aged_70_older + daily_tests_per_100k + life_expectancy + gdp_per_capita_LOG + area_ifr_LOG + completenew_D + Stringency_Treat + Stringency_revised + first.dose.pct + second.dose.pct + I(lag(Stringency_revised, 1:3)) + I(lag(completenew_D, 1:3)) + I(lag(Stringency_Treat, 1:3))  + I(lag(daily_tests_per_100k, 1:3)), data = DF_Full, match.missing = TRUE, verbose = TRUE, qoi = "att", lead = 0:5, outcome.var = "new_cases_UP_100k_change_ihs", placebo.test = T)

cvb_PC_4A <- as.data.frame(get_covariate_balance(Partial_cbps$att,DF_Full, covariates = c("median_age", "gdp_per_capita_LOG", "human_development_index", "v2x_libdem", "v2x_freexp_altinf", "pop.2020_LOG", "aged_70_older", "area_ifr_LOG", "life_expectancy")))

cvb_PC_4B <- as.data.frame(get_covariate_balance(Partial_cbps$att,DF_Full, covariates = c("new_cases_UP_100k_change_ihs", "daily_tests_per_100k", "Stringency_revised", "Stringency_Treat", "completenew_D")))


cvb_PC_1A <-apply(as.data.frame(apply(cvb_PC_1A, 1, abs)),2, mean)
cvb_PC_1B <-apply(as.data.frame(apply(cvb_PC_1B, 1, abs)),2, mean)
cvb_PC_2A <-apply(as.data.frame(apply(cvb_PC_2A, 1, abs)),2, mean)
cvb_PC_2B <-apply(as.data.frame(apply(cvb_PC_2B, 1, abs)),2, mean)
cvb_PC_3A <-apply(as.data.frame(apply(cvb_PC_3A, 1, abs)),2, mean)
cvb_PC_3B <-apply(as.data.frame(apply(cvb_PC_3B, 1, abs)),2, mean)
cvb_PC_4A <-apply(as.data.frame(apply(cvb_PC_4A, 1, abs)),2, mean)
cvb_PC_4B <-apply(as.data.frame(apply(cvb_PC_4B, 1, abs)),2, mean)


cvb_plot <- as.data.frame(rbind(cvb_PC_1A, cvb_PC_1B, cvb_PC_2A, cvb_PC_2B, cvb_PC_3A, cvb_PC_3B, cvb_PC_4A, cvb_PC_4B))

cvb_plot$Refinement <- factor(c("None","None","Mahalanobis","Mahalanobis","PS","PS","CBPS","CBPS"), levels= c("None", "Mahalanobis", "PS", "CBPS"))

cvb_plot$`Variable Type` <- as.factor(c("Fixed","Time-Variant","Fixed","Time-Variant","Fixed","Time-Variant","Fixed","Time-Variant"))


cvb_plot <- melt(cvb_plot, id.vars=c("Refinement", "Variable Type"))
colnames(cvb_plot)[3] <- "Time Period"
colnames(cvb_plot)[4] <- "Average of Standardized Mean Differences"

ggplot(cvb_plot, aes(`Time Period`, `Average of Standardized Mean Differences`, group = interaction(Refinement, `Variable Type`),  color = Refinement, linetype = `Variable Type`)) + geom_hline(yintercept=0.20, color = "black", size=0.9)+ geom_line(size=1.8) + geom_point(size=6.5, shape=18) + theme_bw(base_size = 22) + scale_y_continuous(breaks=seq(0.0, 0.6, 0.1), limits=c(0, 0.6))










